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§ ■ ABSTRACT 
(N 

^ ■ HD 38529 is a post-main sequence G8III/IV star (3.5 Gyr old) with a plan- 

etary system consisting of at least two planets having Msini of 0.8 Mj^p and 
^ . 12.2 Mjup, semimajor axes of 0.13 AU and 3.74 AU, and eccentricities of 0.25 

and 0.35, respectively. Spitzer observations show that HD 38529 has an excess 
^ . emission above the stellar photosphere, with a signal-to-noise ratio (S/N) at 70 

/im of 4.7, a small excess at 33 fim (S/N=2.6) and no excess <30 fim. We 
discuss the distribution of the potential dust-producing planetesimals from the 
. study of the dynamical perturbations of the two known planets, considering in 

particular the effect of secular resonances. We identify three dynamically sta- 
ble niches at 0.4-0.8 AU, 20-50 AU and beyond 60 AU. We model the spectral 
^ , energy distribution of HD 38529 to find out which of these niches show signs of 

\0 . harboring dust-producing plantesimals. The secular analysis, together with the 

SED modeling resuls, suggest that the planetesimals responsible for most of the 
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\^ , dust emission are likely located within 20-50 AU, a configuration that resembles 

O I that of the Jovian planets + Kuiper Belt in our Solar System. Finally, we place 

O ' upper limits (8x10^^ lunar masses of 10 /im particles) to the amount of dust 

that could be located in the dynamically stable region that exists between the 
: two planets (0.25-0.75 AU). 
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Subject headings: circumstellar matter — Kuiper Belt — infrared: stars — plan- 
etary systems — stars: HD 38529 



1. Introduction 



HD 38529 is a pos t-mai n sequence star (G8 III/I V) with an estimated age of 3.5±1 
Gyr (Valenti & Fischer l2005l and Gonzalez et al. 1200 ll ). a distance of 42±2 pc (Ferryman 
et al. Il997). and Tg// = 5697 K, L, = 6.31 Lq, M, = 1.47 Mq [Fe/H] = 0.445 (Valenti 
& Fischer 120051 ). High precision radial velocity monitoring of HD 3 8529 h as led to the 
discovery of a close-in Jupiter-mass planet (HD 38529b - Fischer et a l. 120011) and a second 
more massive and more distant planet (HD 38529c - Fischer et al. |2003| ). Butler et al. 
( 120061 ) report the current estimates of the parameters of the two planets: masses (Msinz) of 
0.8 Mjup and 12.2 Mj^p, semimajor axes of 0.13 AU and 3.74 AU, and eccentricities of 0.25 
and 0.35, for HD 38529b and HD 38529c respectively. 



As part of the Spitzer Legacy Program FEPS (Meyer et al. l2006l ). we searched for debris 
disks around the nine stars in the FEPS sample known from radial velocity (RV) studies to 
have one or more massive planets (Moro-Martm et al. 120071 ). HD 38529 was found to be 
the only star in that sub-sample to have an excess emission above the stellar photosphere, 
with a signal-to-noise ratio (S/N) at 70 fim of 4.7, a small excess at 33 fim (S/N = 2.6) and 
no excess at A < 30 fim (Moro-Martm et al. 120071 ). HD 38529 therefore joined the small 
group of stars known to date that have both IR excess and one or more known planetary 
companions. Table 1 summarizes the properties of these systems, showing a wide diversity 
of planetary architectures. Six of these nine sources, HD 33636, HD 50554, HD 52265, HD 
82943, HD 117176 and HD 128311, are similar to HD 38529 in that their Spitzer observations 
also show an excess at 70 /im but no excess at 24 /im, implying that the bulk of the excess 
emiss ion is a rising from cool material (T<100 K) located mainly beyond 10 AU (Beichman 
et al. l2005al ). This means that these stars not only harbor planets (as inferred from their 
radial velocity observations) but also harbor an outer belt of dust-producing planetesimals 
(responsible for their IR excess) and in this regard they resemble the Solar System in its 
Jovian planets + Kuiper Belt configuration. The other two sources in Table 1 (HD 69830 
and e-Eri) have warm dust and are therefore less relevant to the present discussion. 

In this paper we constrain the distribution of the potential dust-producing planetesimals 
from the study of the dynamical perturbations of the two known planets, considering in 
particular the effect of secular resonances. This allows us to identify three regions where 
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planetesimals could be dynamically stable (§2). After examining the lifetimes of the dust 
grains and concluding that the dust traces the location of the parent planetesimals (§3), we 
show how we can further constrain the location of the dust (and the parent planetesimals) 
by modeling HD 38529 spectral energy distribution (SED - §4). In §5 we discuss together 
the results from the dynamical and SED analysis, and finally §6 summarizes our results. 



2. Possible Location of the Dust-Producing Planetesimals: Effect of 
Gravitational Perturbations by the Planets 

We can identify the possible location of the dust-producing planetesimals by studying 
the effect of the planetary perturbations on the stability of the planetesimals' orbits. Consider 
the orbital parameter space of semimajor axis and eccentricity, (a, e). First, we can eliminate 
planetesimals with orbits that would cross the orbits of the planets, i.e., orbits that have 
either periastron or apoastron within the periastron-apoastron range of each of the planets. 
Second, we note that planetesimals in initially circular orbits would be strongly unstable in 
the vicinity of the orbits of each of the t wo kno wn planets, in a range of semimajor axis, 
Aa ~ ±1.5(mp|g^j^g|-/m^)^/'' (Duncan et al. Il989l ). These two considerations identify several 



regions where planetesimals could not be stable, shown as the grey and red shaded zones 
in Fig. 1. In addition to these perturbations that operate over short timescales, we need to 
consider the effect of perturbations operating over much longer timescales. These are the 
secular perturbations, and as we will see below, they can cause a strong eccentricity excitation 
of the planetesimals, particularly at secular resonance locations, that can significantly shorten 
their lifetime. 

We now calculate the effect of secular perturbations of the two planets on the planetes- 
imals (taken as test particles on initially circular orbits). For the two planets, we assume 
co-pla nar orbits and minimum masses, taking their orbital parameters from Butler et al. 



( 120061 ). For a two planet system, there are two linear modes that excite eccentricities of 
test parti cles; w e follow the Laplace-Lagrange secular perturbation analysis (see Murray and 
Dermott Il999h . The frequencies, gi, phases. Pi, and amplitudes, E^^\ of the two secular 
modes for the HD 38529 system are: 

gi = 0.106arcsec/yr, E^^^ = (6.697 x 10-^ 0.3300), /3i = 0.2271 (1) 

c/2 = 19.7arcsec/yr, E^^^ = (0.2787, -6.411 x 10'^), P2 = 1-594 (2) 

The secular variations of the eccentricity vectors of planets b,c are given by 

Cb cos zub = e[^^ cos{git + /3i) + e[^^ cos{g2t + (32) (3) 
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Cb sin = e[^^ sm{git + /3i) + Ef^ sin(^2t + ^2) (4) 
Cc COS Wc = E2^ cos((7it + (3i) + £'2^^ cos{g2t + (32) (5) 
Cc sin Wc = E2^^ sin((7it + (3i) + £'2^'' sin((72t + (32) (6) 

The fast mode dominates in planet b's eccentricity vector, whereas the slow mode dom- 
inates planet c's eccentricity vector. 

Now consider a test particle (i.e. a planetesimal) in this system. Its secular perturbations 
have two components: a "free" apsidal precession due to the overall quadrupole potential 
(arising from the combined effects of the post-Newtonian stellar gravitational field, and the 
gravitational effects of planetary masses approximated by uniform density circular rings of 
radius equal to their semimajor axes), and a "forced" eccentricity-apsidal perturbation due 
to the planetary secular modes (i.e., due to the secular variations of the elliptical planetary 
orbits). The amplitude of the latter perturbation is a function of the semimajor axis of the 
test particle; at some values of semimajor axis, where the free precession frequency is close 
to a planetary secular frequency, the test particle is subjected to a large amplitude resonant 
eccentricity excitation. To determine the maximum forced eccentricity obtained by a test 



particle, we use the secular resonance analysis given in Malhotra (119981 ) which includes the 
effects of non-linear saturation of the amplitude of the resonantly forced eccentricity from 
each secular mode. The results are shown in Fig. 1: a test particle initially in a circular 
orbit would have its eccentricity excited by the slow mode to the values indicated by the 
blue curve; the green curve is for the eccentricity excited by the fast mode. 

Clearly, very significant secular eccentricity excitation occurs over a wide zone that 
extends to distances much larger than the unstable zones identified above. Remarkably, 
even though the strongly unstable zone of the outer planet extends outward to only about 
5.5 AU, the eccentricity excitation exceeds 0.3 to more than 10 AU. In general, we see that 
the forced eccentricities exceed 0.1 everywhere up to ~ 57 AU. A strong secular resonance 
with the slow mode occurs at a semimajor axis value of ~ 55 AU, where the particles become 
planet-crossing or can even collide with the star; from the secular analysis (Malhotra 1998), 
we calculate that the timescale to excite the eccentricity from ~ to ~ 1 is ~ 1.9 x 10^ years. 
(An analogous phenomenon occurs at the inner edge of the asteroid belt in our Solar System, 
at the location of the so-called z/g secular resonance where asteroidal apsidal precession rates 
are resonant with the 6th secular mode of the Solar System.) Thus, tentatively, we can 
identify three regions that could harbor planetesimal populations (or perhaps small planets) 
in low eccentricity orbits: an inner region 0.4-0.8 AU in-between the two planets' orbits, and 
two outer regions, 20-50 AU and beyond 60 AU. This conclusion is confirmed by numerical 
integrations of test particles as demonstrated in Figures 2 and 3. 
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What does the high eccentricity excitation of test particles imply for dust-producing 
planetesimal populations? Because these are forced eccentricities, neighboring particle orbits 
could have a high degree of apsidal alignment, so that the relative velocities of planetesimals 
would be determined by Kepler shear, and therefore would be small. However, it is likely 
that, in the formation and evolution process, free eccentricities of planetesimal populations 
become comparable to the forced eccentricities, as observed in the solar system's asteroid 
belt and Kuiper Belt. In that case, the forced eccentricities would indeed be a measure of the 
relative velocities for collisions amongst planetesimals in these regions. Thus, if the zones 
identified above are populated by planetesimals, the planetary perturbations would cause 
mutual planetesimal collisions that would result in the production of dust. 

Alternatively, because the primordial disk of HD 38529 was massive enough to form at 
least two giant planets, arguably, these zones could harbor one or more large planetesimals, 
or even sub- Jovian mass planets so far undetected: sufficiently massive planets would tend 
to suppress the secularly forced eccentricities, and clear out planetesimal populations from 
their vicinities, thereby removing sources of dust; while large 1000 km-size planetesimals 
could have stirred up and ground away the planetesimal disk at early times. 



3. A Collision-Dominated System: the Dust Traces the Planetesimals 

To evaluate whether the HD 38529 debris disk is radiation-dominated or collision- 
dominated we consider the lifetimes of the dust grains: 



The collisional lifetime at a distance R from the star can be estimated as t 



col 



Porb/8S(Tc,eo (Backmau & Paresce Il993l ). where Port is the orbital period, S is the 
surface density of the grains and ag^o is the grain geometric cross section. Scr^eo 
is the dimensionless fractional surface area of the disk, of the order of l^dust/'^* ~ 
10-5(5600/T,)3(F7o,d„,t/F7o,*) (Backman & Paresce |l993|) . For HD 38529, F70,, = 17 4 
mJy, F7o,d„st = 66.9 mJy, T, = 5697 K and M, = 1.47 M© (Moro-Martm et al. l2007h . 
so that hdust/'^* = 3.6 xlO~^ and 

t,oi = 3.4 X 10=^P,,, = 2.8 X 10=^( A)3/2^^. (7) 

The Poynting-Robertson (P-R) lifetime at distance R, i.e. the time it takes for the 
dust grain to migrate from R to the star, is given by 

tpR = ^-T- = 710(— )(^— ^)(777) (7^) . _^ , yr (8) 
3 L^, jjm g/cm-^ AU L^, 1 + albedo 
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(Bums, Lamy & Soter, Il979l : Backman & Paresce Il993l ). where b and p are the grain 
radius and density, respectively. More relevant is the time it takes for the particle to 
drift under P-R drag from a populated to a relatively unpopulated region, i.e. the time 



it would take to fill out the gap under P-R drag (tp^). If the scale over which the dust 
density significantly decreases is x% of R, this time is 



t 



fill 

PR 



:i-a 



X 



100' 



tpR = 2.0 X 10^ 



, R , 



(9) 



where we have assumed that x = 10%, L* = 6.3 Lq, p = 2.5 g/cm^, 6 = 10 pm and 
albedo = 0.1. 



If the sources of dust are outside the orbit of the planet, and if P-R drag dominates the 
dynamics, as the dust particles drift inward due to P-R drag they are likely to be scattered out 
of the system when crossing the o rbit of a plan et, cre ating a dust depleted regi on inside th e 
planet's orbit (Roques et al. llDoi Liou & Zook|l999|; Moro-Martin & Malhotra |20oi, IsQOsI) . 
The ejection is very efficient for planets in circular orbits: planets with masses of 3-10 Mjup 
located between 1-30 AU eject >90% of the particles that go past their orbits, while a 1 
Mjup planet at 30 AU ejects > 80% of the particles, and about 50%-90% if located at 1 



AU (Moro-Martin & Malhotra l2005l ): these results are valid for dust particle sizes in the 
range 0.7-135 pm (assuming astronomical silicate composition), and when the central star 
is solar- type. 



On the other hand, if the system is collision-dominated (Krivov et al. l2000l : Dominik & 
Decin l2003l : Wyatt l2005l ). mutual collisions can fragment the larger particles to smaller and 
smaller sizes, until they are blown out from the system by radiation pressure. This means 
that the dust particles would be destroyed before they migrate far from the location of their 
parent bodies under P-R drag, i.e. the dust traces the location of the parent bodies. In this 
case, an inner cavity in the dust density distribution could only arise if the planetesimals 
themselves are confined to a belt, i.e. if there is an inner edge to their spatial distribution. 
This scenario may suggest a massive planet confining the inner edge of the dust-producing 
planetesimals. 

For the HD 38529 system, the estimates in equations [7] and [9] suggest that tcoi << t^l^ 
at all relevant radii, i.e., we are in the collision-dominated regime. This means that the dust 
emission traces the location of the dust-producing planetesimals, with dynamically stable 
niches at 0.4-0.8 AU, 20-50 AU and beyond 60 AU. To find out which of these niches do 
actually show signs of harboring dust-producing plantesimals we turn now to the study of 
the IR excess emission detected by Spitzer. 
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4. SED Modeling 

4.1. Single Temperature Models are Insufficient 

The 70 /iin observations of HD 38529 are not spatially resolved and therefore it is not 
possible to know unambiguously the spatial distribution of the dust. We can learn about the 
characteristic temperature of the dust from the spectral energy distribution (Hillenbrand et 
al. in preparation). From the ratio of the excess emission at 70 fim to that at 33 /im, and 
assuming these wavelengths are at the Wien's tail of the dust blackbody emission, the ratio 
of the excess fluxes yield a characteristic temperature of of 43±4 K. If the blackbody grains 
are in thermal equilibrium, we can calculate the location of the dust from 

Lq Tdust 



(Backman et al. Il993l ): for HD 38529, ^dust = 106±18 AU. The dust mass can then be 
estimated by assuming the dust is in a thin shell at a distance ^dust-, and that the grains are 
spherical with cross section equal to their geometric cross section, so that the total number 
of dust particles is ~(47rR^^^^/7r&^) xr and y^dust ~ 167rRdusf^rfep/3. The optical depth, r, 
can be approximated as hdust/'^* so that, 

6.28 X lQ-^[h^)[^—){^^){^^f (11) 



L^, gem jj,m AU 



(Jura et al. Il995l ). where b and p are the particle's radius and density. It is important to note 
that the above estimate is a lower limit: in the absence of sensitive (sub)millimeter detections, 
no realistic constraints can be made to the dust mass, a significant fraction of which could 
be locked in grains with sizes of ~1 mm that emit efficiently in the (sub)millimeter but 
contribute little to the infrared emission. For HD 38529, Ldust/^* = 3.6x10"^, and using b 
= 10 pm and p = 2.5 g/cm'^ we obtain Mdust > 1.9x10"^ mJ. 

Another estimate of the dust temperature can be made from fitting a photosphere plus 
a single temperature blackbody to the IRS spectrum only; this results in a dust tempera- 
ture of 79 K, which would corresponds to a distance of 31 AU. The discrepancy with the 
characteristic temperature derived from the ratio of the excess emission at 70 pm to that 
at 33 pm indicates that it is not possible to conclude that the dust is confined to a narrow 



^Dus t mass estimates for the KB dust disk range from a total dust mass < 3xlO_Jf M0 (Backman et 
al. llQQSh to - 4x10"" Mq for dust parti cles < 150 (Moro-Martm & Malhotra l2003f ) : with a fractional 
luminosity of Ld„st/L*^10^^-10^^ (Stern 
zodiacal cloud) is estimated to be Ldust/^* 



199< 



The fractional luminosity of the asteroid belt dust (a.k.a 
10~^-10"'^ (Derniott et al. l2Cl02h . 
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ring as implied by the assumption of a single grain temperature. Similarly, Hillenbrand et 
al. (in preparation) find that the SEDs of about 1/3 of the FEPS targets with 70 /im ex- 
cess emission are better fit by multi-temperature rather than single-temperature blackbody 
models. The need for a multi-temperatu re gra in distribution ha s prev iously been found in 
systems hke /3-Pic (e.g. Li & Greenberg Il998l and Telesco et al. l2005l ) and HR4796A (e.g. 
Li & Lunine l2003l ) . and has been unambiguously confirmed by the eight spatially resol ved 
observations in scattered light of debris disks known to date, which led Kalas et al. (120061 ) to 
conclude that debris disks show two basic architectures: 1) belts about 20-30 AU wide and 
with well-defined outer boundaries (HR 4796A, Fomalhaut and HD 139664); and 2) wider 
belts with sensitivity-limited edges implying widths >50 AU (HD 32297, /5-Pic, AU-Mic, HD 
107146 and HD 53143). In the absence of spatially resolved observations, the SED modeling 
of HD 38529 (or any other source under consideration) should enable to explore a range of 
disk architectures. 



4.2. Multi Temperature (SED) Models 

4-2.1. SED Modeling Assumptions 
For the modeling of the SED, we use the radiative transfer code developed by Wolf 



& Hillenbrand (120031 ). Because of the above considerations, we model the dust disk as an 
annulus of inner radius Rj„, outer radius Rout, total dust mass Mdust, and a constant surface 
density (S oc r°, so that the number density, n(r) oc r^^). We assume that the dust grains 



are composed of silicate^ with optical constants from Weingartner & Draine ( 2001 ). For 



the particle sizes we consider two options: 1) a single grain size of 10 /im in radius, and 2) a 
particle size distribution following a power law, n{b) oc where b is the particle radius, q 
= 3.5 (for grains in collisional equilibrium), = 2 fim and bmax = 10 /im. In both cases, 
the radius of 10 was chosen because such a grain radiates efficiently at 70 ^m. Because 
of the last property, this choice of grain size provides a lower limit for Mdust- The value 
selected for bmin corresponds to the "blow-out" size, i.e. the minimum size of the grains that 
can remain bound in the system, given by 

bmin _ Q t;o 2.5g/cm~^ 1 + albedo L^Lq 

(Artymowicz et al. Il989l). For HD 38529, L, = 6.31 Lq, M, = 1.47 Mq and T, = 5697 K; 
using a grain density p = 2.5 g/cm~'^ and albedo = 0.1, we get bmin = 2 /xm . 



^For a study on how the SEDs de pend on the grain composition we refer to Wolf & Hillenbrand ( 2003[ ) 



and Moro-Martm, Wolf & Malhotra (|2005l ). 
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We assume that the dust disk is optically thin with the dust being in thermal equilibrium 
with the stellar radiation field. Only scattering, absorption and reemission of stellar radiation 
by dust grains are taken into account, neglecting scattering and dust heating from the 
infrared radiation produced by the optically thin dust disk. With these assumptions, the 
temperature of the dust grain, for a given dust size and a chemical composition, depends 
only on the distance to the central star §. 

The outer radius of the disk, Rout, cannot be constrained with data currently available. 
Based on scattered-light observations from near by de bris disks, d isk si zes of about 1 00 AU 
to several hu ndred AU are inferre d (Ka las e t al. 



2005: Dent et al. 



2000 



Greaves et al. 



Wilner et al. l2002l : Holland et al. l2003l : Liu l2005l : Metchev et al 
For this reason we consider disk sizes of Rout = 50 AU (Solar System size). Rout 
and Rout = 500 AU. 



2000 



20051 : Ardila et al. \2004 \. 

100 AU 



With the above three values for Rout, assuming a uniform density distribution and with 
the grain size and composition fixed, we then vary Rj„ and M^ust (our only two free param- 
eters) to create a grid of models where we allow Rj„ to vary from the silicate sublimation 
radius (Rsub, where Tgut = 1550 K) to Rout- This accounts for the possibility of having either 
a dust disk of wide radial extent or a narrow ring of dust. 

The observations to be modeled are the ^jsiteer-IR^ synthetic photometric points at 13 
/im, 24 /im and 33 um a nd the Spitzer-MIPS photometric points at 24 fim and 70 fim (from 
Moro-Martin et al. 120071 ). To evaluate whether a particular model is a valid fit to these data 



•^Our SED modeling assumes that the observed dust excess arises from a cicumstellar disk. However, the 
outermost planet has a Msini of 12.2 Mjup, placing it in the boimdary between planets and brown dwarfs. 
Even though the latter could potentially harbor a disk, the prevalence of cold dust in the HD 38529 system 
indicates that this could not be the dominant source of the observed dust because this disk would be located 
near the star (the outermost planet's semimajor axis is 3.74 AU) implying the presence of warm dust that 
is not observed. 

''The IRS spectrum for A < 14.21 fim (short-low module) shows a small offset with respect to the Kurucz 
model. Because the emission at those wavelengths is clearly photospheric, we have multiplied the spectrum 
for A < 14.21 /im by 1.045 to make it coincide with the stellar photosphere (so these wavelengths do not 
dominate the statistics). In addition, the IRS spectrum shows a small discontinuity at 14.21 /j,m, possibly 
because the star was not centered in the slit of the long-low module. To correct for this discontinuity, we have 
multiplied the IRS spectrum for A > 14.21 fim by 1.108. Accounting for this, the corrected 33 /im flux is 95±2 
mjy, where the error is the internal uncertainty. For the calculation of the S /N of the excess we used the total 
uncertainty, obtained from adding in quadrature the internal uncertainty and the calibration uncertainty, the 
latter taken to be 6%. This is an underestimate of the significance of the departure from a pure photosphere 
because this calibration uncertainty is for the overall spectrum, not for individual wavelengths relative to 
each other. 
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points, or if not, to what degree of certainty the model can be excluded, we calculate its 
probability, P(x^ | v\ where v is the number of degrees of freedom. In our case z/ = 2, as 
the only two free parameters are and yidust-, and all the other disk parameters and dust 
properties are fixed to the values given above. The probability is defined so that P(0 | z/) = 
and P(l \ v) = \. Models with P(x^ | ^) > 0.9973 can be excluded with a S-cr certainty, 
while models with P(x^ | v) > 0.683 could be excluded with a l-cr certainty. 

4.2.2. SED Modeling Results 

Figure 4 shows the two-dimensional grids of SED models described above: in green are 
the models with P(x^ | v) < 0.683, i.e. models that are consistent with the observations; 
red represents models with P(x^ | u) > 0.683, i.e. models that could be excluded with a 1-a 
certainty; and black corresponds to models with P(x^ | > 0.9973, i.e. models that can 
be excluded with a S-cr certainty. Some of these SEDs are shown in Fig. 5 together with 
the observations. As can be seen in both figures, the models are degenerate and for a given 
Rout, there are many pairs of Ri„ and Mdust that could fit the observed SED. Because of this 
degeneracy, our main interest lies in identifying the models that can be excluded. These are 
of particular interest because they allow us to identify dust depleted regions whose origin 
can be studied in terms of the overall dynamics of the planetary system (see §2). 

The main results from Fig. 4(a-c) regarding the models that assume grains 10 /im in 
radius are the following 

• For Ro„( = 50 AU, we can exclude models with Rj„ < 5 AU with a certainty of 3-0", 
or we could exclude models with Rin < 14 AU with a certainty of l-cr; the data is 
consistent with models having 15 AU < Rin < 50 AU, i.e., models ranging from a 35 
AU wide disk (extending from 15 AU-50 AU) to a very narrow ring of dust located 
at 50 AU. In other words, the observations require the presence of an inner cavity in 
the dust density distribution of at least 5 AU in radius, with larger cavities consistent 
with the data (the latter implying narrower rings of larger dust mass) . 

• For Rout — 100 AU, none of the Rj„ considered (ranging from the dust sublimation 
radius, Rsub, to 100 AU) can be excluded with a 3-a certainty; while models with Rj„ 
< 7 AU and those with Rin > 40 AU could be excluded to 1-a. Possible fits include 
models with 8 AU < Ri„ < 35 AU. 

• For Rout = 500 AU, again none of the Ri„ considered (from Rgub to 400 AU) can be 
excluded with a 3-a certainty; while models with Ri„ > 14 AU could be excluded to 
1-a. Possible models include those with Rgub ^ Rin ^ 12 AU. 
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We can understand these results as follows. For a constant surface density, S oc r'', 
there is 100 times more dust mass between 50 and 100 AU than there is between 5 and 
10 AU, while the dust temperature is proportional to r~^/^. This results in that there is a 
significant contribution to the 70 /im emission from the largest radii, Kout- To keep the 70 
/im emission constant, the dust surface area or mass at Rout must be similar as Rout varies, 
i.e. TjRout ~ constant, this means that: 1) larger Rout requires a smaller E, and smaller S 
means that we can tolerate smaller inner radii since there is now less dust mass placed there; 
this explains why the SEDs can be fitted with a smaller inner cavity for Rout = 100 AU, and 
with a constant surface density all the way to the sublimation radius in the 500 AU case; 
and 2) smaller Rout requires a larger S, and in order to keep the 24 /im and 33 /im fluxes 
below their upper limits, this requires to increase the inner cavity size. A larger inner cavity 
would also be needed if we had assumed S oc instead of constant. 

In the models shown in Fig. 4(d) we relax the assumption that all the grains are 10 
fim in radius and allow for the presence of smaller grains, with the particle size distribution 
following n{b) oc b'^-^, and with bmin = 2 /im and bmax = 10 /im. We find that most 
of the models can be excluded with a 3-0" certainty, i.e. we can exclude the presence of 
a significant population of small grains inside 100 AU based on the lack of a significant 
continuum emission at A < 30 /im. A similar conclusion has been obtained from several 
other debris disks whose spectroscopy show little or no s olid state features , indic ating that 



the dust grains have sizes > 10 /im (e.g. Jura et al. I2004J : Stapelfeldt et al. 1200J). 



5. Discussion of the Dynamical and the SED Analysis 

5.1. Location of Cold Dust 

The two outermost regions identified in §2 where planetesimals could survive for ex- 
tended periods of time in the presence of the two known radial velocity-detected planets, 
namely 20-50 AU and beyond 60 AU, are in broad agreement with the allowed dust locations 
that result from the modeling of the SED. In particular, we saw that the models with Rout 
= 50 AU predicted an inner cavity of radius > 5 AU; while a range of dust disks with Rout 
= 50 AU and Rj„ > 20 AU could fit the observed SED. From the models with Rout = 100 
AU, we can conclude that it is very unlikely that the observed dust emission arises only 
from planetesimals located beyond 60 AU because the models with Rout = 100 AU and Rin 
= 60 AU, 70 AU, 80 AU and 90 AU have P(x^ | i^) larger than 0.83, 0.87, 0.89 and 0.91, 
respectively A similar conclusion can be drawn from the models with Rout = 500 AU, and Rj„ 
> 60 AU, for which P(x^ | u) > 0.94. Therefore, from the dynamical and SED modeling we 
conclude that the planetesimals responsible for most of the dust emission are likely located 
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within the 20-50 AU region. 

Given that the maximum eccentricities in the 20-50 AU region are moderately high, 
0.25-0.3 (see Figures 1 and 2) - implying coUisional velocities at 30 AU of approximately 
O.Sxvcirc = 2.8 km/s - could a moderate rate of dust production persist in this relatively old 
system? (A "moderate rate" is important in order for the dust production to be long-lived, 
inste ad of being depleted rapidly in the early history of the system). Following Wyatt et 



al. (120071 ). we can estimate the maximum fractional luminosity of the excess, fmax, that 
could originate from a planetesimal belt of a given age that is evolving in quasi-steady state. 
Wyatt's analytical model assumes that the planetesimals and the dust are in coUisional 
equilibrium, and that their size distribution follows a continuous power law of index -3.5. 
For HD 38529, and given the results above, we assume that the planetesimal belt extends 
from 20-50 AU, with a mean planetesimal eccentricity of 0.3. Other parameters in the 
model are the diameter of the largest planetesimal in the cascade, 2000 km, and the specific 
incident energy required to catastrophically destroy a planetesimal, 200 J/kg. At the age of 
the HP 38 529 system, 3.5 Gyr, the model predicts fmax = 2.26x10"^ (see eq. [20] in Wyatt 



et al. 120071 ). This value is smaller than the observed fractional luminosity of the excess, fobs 
= 3.6x10"^ (see §3). However, because there are two orders of magnitude uncertainty in the 
estimate of fmax, we cannot reject a scenario in which the dust observed is the result of the 
steady grinding down of planetesimals. If we were to assume that the above value of fmax is 
correct, an estimate of the timescale over which the fractional luminosity can be mantained 
at the value of fobs is given by tage^ fmax/ fobs = 220 Myr. In this case, a possible scenario 
could be that the stable region beyond 60 AU supplies some planetesimals that drift into 
the 20-50 AU region by non-gravitational effects. 



5.2. Upper Limits on Warm Dust 

In §2 we identified a modestly stable small zone in-between the two known planet s (0.4 - 



0.8 AU), that has also been identified in numerical simulations by Barnes & Raymond (12004 ): 



our secular perturbation analysis provides a theoretical explanation for those numerical sta- 



bility results. Raymond & Barnes (120061 ) consider terrestrial planet accretion in this zone. 
They conclude that HD 38529 is likely to support an asteroid belt and perhaps Mars-sized 
planets, but not larger planets, because the potential feeding zone for the accretion of a ter- 
restrial planet would be limited by the high eccentricities of the planetesimals in this region. 
The lack of an IR excess at wavelengths shorter than 30 /im allows us to place an upper limit 
on the amount of warm (asteroidal) dust that could be located in this inner region. We use 
the IRAC 5/im and 8/im, IRS 13/im and MIPS 24/im photometric measurements, and we 
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assume 10/xm size silicate grains (optical constants taken from Weingartner & Draine |200l|) . 
We find that a 3-cr upper limit to the dust mass in this potential asteroid belt is 3xl0~^^ 
Mq or 10~^ M®. For comparison, the mass estimate for the zo diacal cloud in the terrestrial 
planet region of the Solar System is 3xl0~^° (Hahn et al. 120021 ). i.e. 330 times smaller 
than the estimated upper limit of warm dust in HD 38529. 



6. Conclusions 



HD 38529 harbors a planetary system consisting of at least two planets (with Msin« of 
0.8 Mjup and 12.2 Mj^p, semimajor axes of 0.13 AU and 3.74 AU, and eccentricities of 0.25 
and 0.35) and a likely population of dust-producing planetesimals that are responsible for 
the 70 /zm excess emission detected by Spitzer. Using analytical and numerical dynamical 
analysis, in this paper we have constrained the distribution of the potential dust-producing 
planetesimals from the study of the dynamical perturbations of the two known planets, 
considering in particular the effect of secular resonances. A dust disk inner edge at 5.5 
AU would naturally arise from the gravitational scattering of planetesimals and dust grains 
by the outermost planet. We show that larger inner cavities in the dust disk, that would 
be consistent with the observed SED, can be created due to the secular effects that arise 
from the interaction between the two massive planets. From the analysis of the secular 
perturbations we identify three regions that could harbor planetesimal populations in low 
eccentricity orbits (where the planetesimals could be long-lived): 0.4-0.8 AU, 20-50 AU and 
beyond 60 AU. From the modeling of the SED we conclude that the planetesimals responsible 
for most of the dust emission observed by Spitzer are likely located within the 20-50 AU 
region. In this regard, HD 38529 resembles the configuration of the Solar System's Jovian 
planets + Kuiper Belt (KB) . The SED models give a dust mass estimate of 1-5x10"^° Mq of 
10 yum particles. The presence of a significant population of small grains inside 100 AU (with 
the particle size distribution following n{b) oc b~^'^ , and with bmin = 2 /im and bmax = 10 
/im) is excluded to a 3-0" certainty level based on the lack of a significant continuum emission 
at A < 30 /xm. We do not find evidence of dust emission within the innermost region, with a 
3-0" dust mass upper limit of 10~^ (in 10 fim grains), suggesting any remnant dust belt 
would have a mass smaller than 330 times that in the Solar System's zodiacal cloud. 

The SED models are degenerate. We need to break this degeneracy to get a better 
understanding of this planetary system, in particular, of how the spatial distribution of the 
dust and the planetesimals are affected by the gravitational perturbations of the two planets. 
This requires spatially resolved images to constrain the disk sizes, and/or high resolution 
spectroscopy observations to look for spectral features that could constrain the grain size 
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and composition, and/or accurate photometric points in the 33 //m-70 /im range and in the 
sub- mm to better determine the shape of the SED. 
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wheieYjm^dust and are the dust excess and photospheric flux at 70 /im and T^, is the stellar temperature (Bryden et 
al. (l2006l ). For HD 69830, the fractional luminosity of the excess is calculated by integrating the excess and photospheric 
emission beteween 7 and 35 /xm. References are: (1) Moro-Martm et al. (2007); (2) Beichman et al. (2005a); (3) Beichman 
et al. (2005b) and (4) Greaves et al. (1998, 2005). 
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Fig. 1. — Test particle orbits are planet-crossing (hence unstable) in the red-shaded zone, 
and strongly unstable in the grey-shaded zone (owing to overlapping first order mean motion 
resonances). The secular modes of the two planets, HD38529b and HD38529c, excite the 
eccentricities of circular test particle orbits: the maximum eccentricity due to the fast mode 
is shown by the green curve, while that due to the slow mode is shown by the blue curve. 
The inner region, interior to the outer planet, is shown in the upper panel. The lower panel 
shows that the planetary perturbations are very wide ranging: secular eccentricity excitation 
exceeds 0.1 to nearly 60 AU; the sharp peak at 55 AU is due to a resonance with the slow 
mode. 




Fig. 2. — Numerical simulations of 300 test particles in the HD 38529 system of two planets. 
The simulations include 100 particles uniformly spaced between 0.01-5 AU, and 200 particles 
uniformly spaced between 5-75 AU, all in initially circular orbits coplanar with the planets. 
The angular elements were chosen randomly between and 27r. Particles were removed if they 
approached the star closer than 0.005 AU or approached a planet closer than the Hill radius 
of the pla net. T he orbits were integrated for 200 Myr using a symplectic integrator (Wisdom 
& Holman ll99ll ). Overall, the simulations confirm the results from the secular analysis. They 
differ in that the maximum eccentricity in the 20-50 AU region calculated from the secular 
analysis is smaller than that found in the numerical integrations; this is because the latter 
include secular and non-secular perturbations (e.g. mean motion resonances), and the test 
particles have non-zero initial eccentricities and inclinations. 
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Fig. 3. — Orbital evolution over the 200 Myr integration period of a test particle near the 
secular resonance. The quantity dw is the difference between the longitude of periastron of 
the particle and that of the outermost planet. The timescale to excite the eccentricity is 
about 50 Myr. 



- 23 - 



3 10 - (a) R<,^t = 50 AU 
grain size: lO/xm 

2.5 10-^° - 



2 10 



10 ' - 




Ik 8 10 



- 6 10 



4 10 



2 10 



(b) R„,t = 100 AU 
. grain si/e: lO^ni 



10 20 30 40 50 
Rm (AU) 



20 40 60 80 100 
Rm (AU) 



2 10 



1.5 10 



10 



5 10 



2 10 



(c) R„,t = 500 AU 
grain size: lOJm 



1.5 10 - 



10 



- 5 10 



(d) R„^t = 100 AU 
grain sizes: 2.2— 10/j.m 



;t:::i,m,;,t 



100 200 300 400 500 

Rm (AU) 



20 40 60 80 100 

Rm (AU) 



Fig. 4. — Each point of these two dimensional grids represents a modeled SED, where Rj„ 
and Mdust are the two free parameters. Panels (a), (b) and (c) correspond to models with a 
single grain size of 10 fim and with Kout = 50 AU, 100 AU and 500 AU, respectively. Panel 
(d) corresponds to models with a distribution of grain sizes given by n{b) oc b^^'^ , with bmm 
= 2 fim and b^ax = 10 /xm. The models in green are those with P(x^ | z/) < 0.683; red in 
panels (a)-(c) represents models with P(x^ | t^) > 0.683, i.e. models that can be excluded 
with 1-cr certainty; and black corresponds to models with P(x^ | i^) > 0.9973, i.e. models 
that are excluded with 3-a certainty. The red dots in panel (d) have P(x^ I ^) > 0.988, i.e., 
they are all close to be excluded with a 3-0" certainty. 



Fig. 5. — Observed and modeled SEDs for HD 38529. The dotted line is the Kurucz modeL 
The black thick line is the IRS low-resolution spectrum. The photometric points are iden- 
tified as follows: black circles are the new Spitzer observations (IRAC, MIPS and synthetic 
photometry from IRS); red diamonds are IRAS observations. In all cases, the error bars 
correspond to 1-cr uncertainties. Upper limits are represented by triangles and are given 
when F/AF < 3 and placed at F + 3x AF if F > 0, or 3x AF if F < 0. The colored con- 
tinuous lines in the main panel show three sets of models that fit the observations with a 
probability < 0.68 (corresponding to the models represented in green in Fig. 4). The models 
include the emission from the photosphere and from a dust disk composed of astronomical 
silicate grains 10 fim in radius. We assume the dust disk extends from Rj„ to Rout with a 
constant surface density. We consider three values for Rout- 50 AU (red), 100 AU (blue) and 
500 AU (green). Ri„ and Mdust are allowed to vary. The insert at the lower left shows the 
most relevant excluded models. The solid line represents models excluded with a certainty 
of 3-0", while the dashed line corresponds to 1-cr. The model shown in light blue gives an 
upper limit to the amount of warm dust located between 0.25-0.75 AU. 



